##############################################################################
## File-Name:       MPDS2012b_makeSEandLogitSEdata.R
## Date:            2012-12-10
## Author:          Ken Benoit
## Email:           kbenoit@lse.ac.uk                                      
## Purpose:         Generate Benoit, Laver and Mikhaylov (2009) standard
##                  errors and Lowe et al (2011) scales and SEs for 
##                  the latest (2012b) release of CMP data
## Data Used:       MPDataset_MPDS2012b.dta
##                  (from https://manifestoproject.wzb.eu)
## Packages Used:   foreign, combinat
## Output Files:    MPDS2012b_withSEandLogitSE.Rdata
##                  MPDS2012b_withSEandLogitSE.dta
##                  MPDS2012b_withSEandLogitSE.csv
## Machine:         Ken Benoit's Macbook Air
##############################################################################

rm(list = ls())  # clear

require(foreign)
cmp <- read.dta("MPDataset_MPDS2012b.dta", convert.factors=FALSE)

    
# convert missing data codes to NA
cmp$peruncod[round(cmp$peruncod,2)==99.99] <- NA
cmp$total[round(cmp$total,0)==9999] <- NA
cmp$testresult[round(cmp$testresult,0)==999] <- NA
cmp$pervote[round(cmp$pervote,0)==999] <- NA
cmp$absseat[round(cmp$absseat,0)==999] <- NA
cmp$totseats[round(cmp$totseats,0)==999] <- NA
cmp$presvote[round(cmp$presvote,0)==999] <- NA
cmp$presvote[round(cmp$presvote,0)==998] <- NA

# quick check on check per sums
# NOTE: in the 2012b dataset this indicates that percentages only sum to 100%
#       when the extended categories are also used
cmporig <- cmp
diff <- 100 - apply(cmporig[,21:131],1,sum)
plot(diff)
diff.nonextended <- 100 - apply(cmporig[,c(21:76,131)],1,sum)
plot(diff.nonextended)

# append a list of per categories with SEs
totalIndex <- which(names(cmp)=="total")
perStartIndex <- which(names(cmp)=="per101")
perEndIndex <- which(names(cmp)=="peruncod")
rileIndex <- which(names(cmp)=="rile")
cmp <- data.frame(cmp[,1:(perStartIndex-1)], total=cmp[,"total"], cmp[,perStartIndex:perEndIndex],
                  cmp[,(perEndIndex+2):ncol(cmp)])
cmp <- data.frame(cmp, total=cmp[,"total"], matrix(NA,nrow=nrow(cmp),ncol=ncol(cmp[perStartIndex:perEndIndex])))
names(cmp)[which(names(cmp)=="X1"):ncol(cmp)] <-
  paste(names(cmp)[(perStartIndex:perEndIndex)+1], "SE", sep="")
perStartIndex <- which(names(cmp)=="per101")
perEndIndex <- which(names(cmp)=="peruncod")

### this eliminates 158 cases
# drop if peruncod==NA
# cmp <- cmp[!is.na(cmp$peruncod),]
### this flags 141 Sweden and Norway cases as having their peruncod computed
### maybe not right but hey, then these 141 cases are not right either!
cmp$peruncodcomputed <- is.na(cmp$peruncod)
cmp$peruncod[cmp$peruncodcomputed] <- 100 - rowSums(cmp[cmp$peruncodcomputed,perStartIndex:(perEndIndex-1)])
cmp$peruncod <- apply(as.matrix(cmp$peruncod),1,function(x) max(x,0))
cmp$dropmeNwasNA <- is.na(cmp$total)

## bootstrap the manifesto counts and compute percentages
nman <- nrow(cmp)   # number of manifestos in dataset
nrepl <- 1000       # number of replications required
n <- cmp$total      # take totals from CMP total
p <- cmp[,perStartIndex:perEndIndex]/100 # take given probabilities from CMP per values

# Compute the SE from replicates
# vectorizing and using apply means Look Ma, No Loop
cataggregate <- function(charlist, cmprow) {
  temparray <- cmprow[,charlist[1]] 
  if (length(charlist) > 1) {
    for (i in 2:length(charlist)) {
      temparray <- temparray + cmprow[,charlist[i]]
    }
  }
  return(temparray)
}
require(combinat)
cmpSE <- function(cmprow, replicates=100,
                  scale=c("none", "saliency", "logitPos", "logitImp"),
                  charlist.L=NULL, charlist.R=NULL, offset=0.5) {
  scale.type <- match.arg(scale)
  p <- matrix(cmprow[-1]/100, nrow=1)  # the per categories
  n <- as.numeric(cmprow[1])           # total
  if(is.na(n)) n <- 200                # fix to something closed to the overall mean if n==NA to stop it breaking
  # if this is one of those "all missing" rows, return vector of NAs for per categories if "none" 
  # or just a single NA if one of the specific variables and one of the Lowe et al types
  if (sum(is.na(cmprow))==length(cmprow)) if (scale.type=="none") return(rep(NA, length(cmprow)-1)) else return(NA) 
  tempdraw <- rmultinomial(rep(n,replicates), p)/n*100
  colnames(tempdraw) <- names(cmprow)[-1]
  if (scale.type=="none") { 
    result <- apply(tempdraw, 2, sd)
  } else {
    if (is.null(charlist.L) & is.null(charlist.R)) {
      stop("Must specify at least L and R; make charlist.L=NA if only added per categories like markeco")
    }
    if (is.na(charlist.L[1]))  L<-0 else L <- cataggregate(charlist.L, tempdraw)
    R <- cataggregate(charlist.R, tempdraw)
    if (scale.type=="saliency")  {
      result <- sd(R-L)
    }
    else {
      offsetrescaled <- offset / n * 100  # because we are using %s not Ns
      if (scale.type=="logitPos") {
        result <- sd(log((R+offsetrescaled)/(L+offsetrescaled)))
      } else if (scale.type=="logitImp") {
        result <- sd(log((R+L+(offsetrescaled*2))/n))
      } else stop(paste(scale.type, ": Illegal value of scale", sep=""))
    }
  }
  result
}


# simulate per category SEs
nrep <- nrepl
perseStartIndex <- which(names(cmp)=="per101SE")
perseEndIndex <- which(names(cmp)=="peruncodSE")
cmp[,perseStartIndex:perseEndIndex] <-
  t(apply(cmp[,(perStartIndex-1):perEndIndex], 1, cmpSE, replicates=nrep))

# original rile SE
r <- c("per104","per201","per203","per305","per401","per402","per407","per414",
       "per505","per601","per603","per605","per606")
l <- c("per103","per105","per106","per107","per202","per403","per404","per406",
       "per412", "per413","per504","per506","per701")
cmp$rileSE <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                    cmpSE, scale="saliency", replicates=nrep,
                    charlist.L=l, charlist.R=r)
# original planeco SE
l <- NA
r <- c("per403","per404","per412")
cmp$planecoSE <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                       cmpSE, scale="saliency", replicates=nrep,
                       charlist.L=l, charlist.R=r)
# original markeco SE
r <- c("per401","per414")
l <- NA
cmp$markecoSE <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                       cmpSE, scale="saliency", replicates=nrep,
                       charlist.L=l, charlist.R=r)
# original welfare SE
r <- c("per503", "per504")
l <- NA
cmp$welfareSE <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                       cmpSE, scale="saliency", replicates=nrep,
                       charlist.L=l, charlist.R=r)
# original europ SE
r <- "per108"
l <- "per110"
cmp$europSE <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                       cmpSE, scale="saliency", replicates=nrep,
                       charlist.L=l, charlist.R=r)


cmp[cmp$dropmeNwasNA,(perseStartIndex:perseEndIndex)] <- NA
save(cmp, file="cmpwithSE.Rdata")


#### Now compute Lowe et al 2011 logit scales and SEs

ci.lo <- function(x) { quantile(x, probs=.025) }
ci.hi <- function(x) { quantile(x, probs=.975) }
make.pos <- function(cmpdata, charlist.R, charlist.L, offset=.5) {
  R <- cataggregate(charlist.R, cmpdata)
  L <- cataggregate(charlist.L, cmpdata)
  log((R+offset)/(L+offset))
}
make.imp <- function(cmpdata, charlist.R, charlist.L, n=100, offset=.5) {
# because N is default 100, assumes R and L are per
  R <- cataggregate(charlist.R, cmpdata)
  L <- cataggregate(charlist.L, cmpdata)
  log((R+L+(offset*2))/n)
}

rescaledoffset <- 0.5 / cmp$total * 100
# overwrite cases where total is missing with total=200 (7 of these)  
rescaledoffset[is.na(cmp$total)] <- 0.5 / 200 * 100

####
#### SOME TESTS TO MAKE SURE THE % instead of N is working OK with rescaled offset
####
r <- "per302"
l <- "per301"
cmp$decentralization        <- make.pos(cmp, r, l, offset=rescaledoffset)
cmp$decentralization.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
## Position
(N <- cmp$total[814])
(off <- 0.5 / N * 100)
50/81
(R <- cmp$per302[814])
(L <- cmp$per301[814])
log((R+off)/(L+off))
log(0+50/81) - log(3.7+50/81)
log(0+.5) - log(round(3.7/100*81)+.5)
cmp$decentralization[814]
cmp$decentralization.SE[814]
## Importance
log((R+2*off+L)/100)
log(R + L + 100/N) - log(100)
cmp$decentralization.imp[814]
cmp$decentralization.impSE[814]


r <- "per102"
l <- "per101"
cmp$foreignalliances        <- make.pos(cmp, r, l, offset=rescaledoffset)
cmp$foreignalliances.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
#cmp$foreignalliances.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$foreignalliances.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$foreignalliances.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$foreignalliances.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per104"
l <- "per105"
cmp$militarism        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$militarism.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$militarism.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$militarism.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                               cmpSE, scale="logitPos", replicates=nrepl,
                               charlist.L=l, charlist.R=r)
cmp$militarism.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$militarism.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                               cmpSE, scale="logitImp", replicates=nrepl,
                               charlist.L=l, charlist.R=r)
r <- "per109"
l <- "per107"
cmp$internationalism        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$internationalism.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$internationalism.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$internationalism.SE     <-apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$internationalism.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$internationalism.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per110"
l <- "per108"
cmp$logeu        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$logeu.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$logeu.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$logeu.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$logeu.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$logeu.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per204"
l <- "per203"
cmp$constitutionalism        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$constitutionalism.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$constitutionalism.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$constitutionalism.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$constitutionalism.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$constitutionalism.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per302"
l <- "per301"
cmp$decentralization        <- make.pos(cmp, r, l, offset=rescaledoffset)
cmp$decentralization.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
#cmp$decentralization.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$decentralization.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$decentralization.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$decentralization.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per407"
l <- "per406"
cmp$protectionism        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$protectionism.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$protectionism.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$protectionism.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$protectionism.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$protectionism.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per414"
l <- "per409"
cmp$keynesian        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$keynesian.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$keynesian.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$keynesian.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$keynesian.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$keynesian.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per601"
l <- "per602"
cmp$nationalism        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$nationalism.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$nationalism.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$nationalism.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$nationalism.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$nationalism.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per603"
l <- "per604"
cmp$tradmorality        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$tradmorality.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$tradmorality.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$tradmorality.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$tradmorality.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$tradmorality.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per608"
l <- "per607"
cmp$multicult        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$multicult.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$multicult.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$multicult.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$multicult.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$multicult.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per702"
l <- "per701"
cmp$laborpolicy        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$laborpolicy.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$laborpolicy.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$laborpolicy.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$laborpolicy.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$laborpolicy.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- c("per403","per412","per413","per415")
l <- c("per401","per402")
cmp$freemarket        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$freemarket.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$freemarket.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$freemarket.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$freemarket.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$freemarket.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- c("per410")
l <- c("per501","per416")
cmp$environment        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$environment.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$environment.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$environment.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$environment.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$environment.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- c("per401","per402","per407","per414","per505")
l <- c("per403","per404","per406","per412","per413","per504","per506","per701")
cmp$stateconomy        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$stateconomy.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$stateconomy.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$stateconomy.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$stateconomy.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$stateconomy.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- c("per505","per507")
l <- c("per504","per506")
cmp$stateservices        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$stateservices.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$stateservices.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$stateservices.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$stateservices.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$stateservices.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- c("per104","per201","per203","per305","per401","per402","per407","per414",
       "per505","per601","per603","per605","per606")
l <- c("per103","per105","per106","per107","per202","per403","per404","per406",
       "per412", "per413","per504","per506","per701")
cmp$logrile        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$logrile.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$logrile.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$logrile.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$logrile.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$logrile.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- c("per401","per414")
l <- c("per403","per404","per412")
cmp$logplaneco        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$logplaneco.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$logplaneco.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$logplaneco.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$logplaneco.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$logplaneco.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per505"
l <- "per504"
cmp$logwelfare        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$logwelfare.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$logwelfare.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$logwelfare.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$logwelfare.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$logwelfare.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- c("per104","per201","per203","per305","per601","per603","per605","per606")
l <- c("per103","per105","per106","per107","per202")
cmp$loglibcons        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$loglibcons.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$loglibcons.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$loglibcons.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$loglibcons.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$loglibcons.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)

r <- "per507"   #education expenditure limitation
l <- "per506"   #education expenditure expansion
cmp$education        <- make.pos(cmp, r, l, offset=rescaledoffset)
#cmp$education.lo95   <- make.pos(cmp, r, l, ci.lo)
#cmp$education.hi95   <- make.pos(cmp, r, l, ci.hi)
cmp$education.SE     <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitPos", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)
cmp$education.imp    <- make.imp(cmp, r, l, offset=rescaledoffset)
cmp$education.impSE  <- apply(cmp[,(perStartIndex-1):perEndIndex], 1,
                                     cmpSE, scale="logitImp", replicates=nrepl,
                                     charlist.L=l, charlist.R=r)


save(cmp, file="MPDS2012b_withSEandLogitSE.Rdata")

# overwrite logit scale SEs with NAs for those 7 manifestos where total=NA
# SEindex <- grep("SE", names(cmp))
# cmp[cmp$dropmeNwasNA,SEindex] <- NA
# drop columns total.1 and dropmeNwasNA
if (length(grep("total.1", names(cmp)))>0) { cmp <- cmp[,-grep("total.1", names(cmp))] }
if (length(grep("dropmeNwasNA", names(cmp)))>0) { cmp <- cmp[,-grep("dropmeNwasNA", names(cmp))] }
if (length(grep("peruncodcomputed", names(cmp)))>0) { cmp <- cmp[,-grep("peruncodcomputed", names(cmp))] }

write.dta(cmp, file="MPDS2012b_withSEandLogitSE.dta", version=10)
write.csv(cmp, "MPDS2012b_withSEandLogitSE.csv")

